From XML data#
Example on how to load your data if it is currently in xml format
import json
import dolfin
from fenics_plotly import plot
import pulse
try:
from dolfin_adjoint import Constant, DirichletBC, Function, Mesh, interpolate
except ImportError:
from dolfin import Mesh, Function, Constant, DirichletBC, interpolate
# Mesh
mesh = Mesh(pulse.utils.mpi_comm_world(), "data/mesh.xml")
# Marker functions
facet_function = dolfin.MeshFunction("size_t", mesh, "data/facet_function.xml")
marker_functions = pulse.MarkerFunctions(ffun=facet_function)
# Markers
with open("data/markers.json", "r") as f:
markers = json.load(f)
# Fiber
fiber_element = dolfin.VectorElement(
family="Quadrature",
cell=mesh.ufl_cell(),
degree=4,
quad_scheme="default",
)
fiber_space = dolfin.FunctionSpace(mesh, fiber_element)
fiber = Function(fiber_space, "data/fiber.xml")
microstructure = pulse.Microstructure(f0=fiber)
# Create the geometry
geometry = pulse.HeartGeometry(
mesh=mesh,
markers=markers,
marker_functions=marker_functions,
microstructure=microstructure,
)
activation = Function(dolfin.FunctionSpace(geometry.mesh, "R", 0))
activation.assign(Constant(0.0))
matparams = pulse.HolzapfelOgden.default_parameters()
material = pulse.HolzapfelOgden(
activation=activation,
parameters=matparams,
f0=geometry.f0,
)
# LV Pressure
lvp = Constant(0.0)
lv_marker = markers["ENDO"][0]
lv_pressure = pulse.NeumannBC(traction=lvp, marker=lv_marker, name="lv")
neumann_bc = [lv_pressure]
# Add spring term at the base with stiffness 1.0 kPa/cm^2
base_spring = 1.0
robin_bc = [pulse.RobinBC(value=Constant(base_spring), marker=markers["BASE"][0])]
# Fix the basal plane in the longitudinal direction
# 0 in V.sub(0) refers to x-direction, which is the longitudinal direction
def fix_basal_plane(W):
V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
bc = DirichletBC(V.sub(0), Constant(0.0), geometry.ffun, markers["BASE"][0])
return bc
dirichlet_bc = [fix_basal_plane]
# You can also use a built in function for this
# from functools import partial
# dirichlet_bc = partial(pulse.mechanicsproblem.dirichlet_fix_base_directional,
# ffun=geometry.ffun,
# marker=geometry.markers["BASE"][0])
# Collect boundary conditions
bcs = pulse.BoundaryConditions(
dirichlet=dirichlet_bc,
neumann=neumann_bc,
robin=robin_bc,
)
# Create the problem
problem = pulse.MechanicsProblem(geometry, material, bcs)
problem.solve()
2024-06-27 20:58:22,336 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:58:22,337 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:58:22,338 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:58:22,338 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:58:22,339 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:58:22,339 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:58:22,340 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:58:22,340 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:58:22,341 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:58:22,342 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:58:22,342 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:58:22,343 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:58:22,343 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:58:23,110 [734] DEBUG UFL_LEGACY: Blocks of each mode:
99 full
2024-06-27 20:58:23,185 [734] DEBUG UFL_LEGACY: Blocks of each mode:
2024-06-27 20:58:23,393 [734] DEBUG UFL_LEGACY: Blocks of each mode:
18 full
2024-06-27 20:59:13,700 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:59:13,701 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:59:13,702 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:59:13,702 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:59:13,703 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:59:13,703 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:59:13,704 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:59:13,704 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:59:13,705 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:59:13,705 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:59:13,706 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:59:13,708 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:59:13,708 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:59:13,848 [734] DEBUG UFL_LEGACY: Blocks of each mode:
10 full
2024-06-27 20:59:13,923 [734] DEBUG UFL_LEGACY: Blocks of each mode:
3 full
2024-06-27 20:59:14,156 [734] DEBUG UFL_LEGACY: Blocks of each mode:
3 full
(2, True)
# Solve the problem
pulse.iterate.iterate(problem, (lvp, activation), (1.0, 0.2))
2024-06-27 20:59:15,291 [734] INFO pulse.iterate: Iterating to target control (f_26)...
2024-06-27 20:59:15,292 [734] INFO pulse.iterate: Current control: f_26 = 0.000,0.000
2024-06-27 20:59:15,292 [734] INFO pulse.iterate: Target: 1.000,0.200
*** Warning: PETSc SNES solver diverged in 1 iterations with divergence reason DIVERGED_DTOL.
*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.
*** Warning: PETSc SNES solver diverged in 3 iterations with divergence reason DIVERGED_DTOL.
([Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 76),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 259),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 316),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 365),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 414),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 455),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 496),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 545),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 618),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 676)],
[(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 72),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 74)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 255),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 257)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 312),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 314)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 361),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 363)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 410),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 412)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 451),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 453)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 492),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 494)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 541),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 543)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 614),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 616)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 672),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 674))])
# Get the solution
u, p = problem.state.split(deepcopy=True)
volume = geometry.cavity_volume(u=u)
print(f"Cavity volume = {volume}")
2024-06-27 20:59:21,198 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:59:21,199 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:59:21,200 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:59:21,200 [734] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-06-27 20:59:21,336 [734] DEBUG UFL_LEGACY: Blocks of each mode:
1 full
Cavity volume = 1.1189963024197862
# Move mesh accoring to displacement
u_int = interpolate(u, dolfin.VectorFunctionSpace(geometry.mesh, "CG", 1))
mesh = Mesh(geometry.mesh)
dolfin.ALE.move(mesh, u_int)
fig = plot(geometry.mesh, color="red", show=False)
fig.add_plot(plot(mesh, show=False))
fig.show()